function SAT_SST_Step_03_get_scaling_validation(P)

    Out      = read_data(P);
    
    % *********************************************************************
    % Compute output arguments
    % *********************************************************************
    clear('R','E','Validation_data','Parameters')
    for ct_x = 1:1:36
        for ct_y = 1:1:18
            
            % Out_grid   = prepare_grid_data(Out,ct_x,ct_y,P);
            clear('Out_training','Out_validation')
            [Out_training,Out_validation] = prepare_grid_data(Out,ct_x,ct_y,P);
            
            if ~isempty(Out_training)
                
                disp('Get data finished')
                
                for ct_val = 1:numel(Out_training)
                
                    disp(num2str([ct_val numel(Out_training)],'%3.0f / %3.0f'))
                    
                    clear('T_a_training','dT_a_dt_training')
                    T_a_training     = BB98d_forcing_month2high_reso(Out_training{ct_val}.AT_itp);
                    T_a_training     = T_a_training - nanmean(T_a_training(1:100));
                    dT_a_dt_training = [(T_a_training(2:end) - T_a_training(1:end-1)) ./ (86400/2); 0];

                    % Fit parameters assuming no seasonality --------------
                    % x_fit = BB98d_fit_parameters(T_a,dT_a_dt,T_o,l_use_loss,do_season,do_spec)
                    x_fit = BB98d_fit_parameters(T_a_training,dT_a_dt_training,...
                        Out_training{ct_val}.SST,Out_training{ct_val}.l_use_loss,P);
                    Parameters{ct_x,ct_y}.parameters(:,ct_val) = x_fit;

                    clear('T_a_validation','dT_a_dt_validation')
                    T_a_validation     = BB98d_forcing_month2high_reso(Out_validation{ct_val}.AT_itp);
                    T_a_validation     = T_a_validation - nanmean(T_a_validation(1:100));
                    dT_a_dt_validation = [(T_a_validation(2:end) - T_a_validation(1:end-1)) ./ (86400/2); 0];
                
                    P_fit = BB98d_predict(T_a_validation,dT_a_dt_validation,x_fit);
            
                    Out_validation{ct_val}.AT_cowtan = Out_validation{ct_val}.AT * (0.86 - 0.25*Out_validation{ct_val}.lf/100);
                    Out_validation{ct_val}.AT_scl    = P_fit.T_o_m;
                    Out_validation{ct_val}.AT_scl(isnan(Out_validation{ct_val}.AT)) = nan;

                    Validation_data{ct_x,ct_y}.training   = Out_training;
                    Validation_data{ct_x,ct_y}.validation = Out_validation;
                    
                    % Calculate R -----------------------------------------
                    clear('at','at_scl','at_cow','sst')
                    sst    = Out_validation{ct_val}.SST       + Out_validation{ct_val}.AT  * 0;
                    at     = Out_validation{ct_val}.AT        + Out_validation{ct_val}.SST * 0;
                    at_cow = Out_validation{ct_val}.AT_cowtan + Out_validation{ct_val}.SST * 0;
                    at_scl = Out_validation{ct_val}.AT_scl    + Out_validation{ct_val}.SST * 0;
                    
                    if P.do_season == 1
                        intv = 12;
                    else
                        intv = 1;
                    end
                    
                    R{ct_x,ct_y}.R_raw(ct_val) = CDC_corr(CDC_demean(sst,1,intv),CDC_demean(at,1,intv));
                    R{ct_x,ct_y}.R_scl(ct_val) = CDC_corr(CDC_demean(sst,1,intv),CDC_demean(at_scl,1,intv));
                    
                    % Calculate Error -------------------------------------
                    E{ct_x,ct_y}.E_raw(ct_val) = nanmean((CDC_demean(sst,1,intv) - CDC_demean(at,1,intv)).^2);
                    E{ct_x,ct_y}.E_scl(ct_val) = nanmean((CDC_demean(sst,1,intv) - CDC_demean(at_scl,1,intv)).^2);
                    E{ct_x,ct_y}.E_cow(ct_val) = nanmean((CDC_demean(sst,1,intv) - CDC_demean(at_cow,1,intv)).^2);
                end  
                
            else
                R{ct_x,ct_y} = [];
                E{ct_x,ct_y} = [];
                Validation_data{ct_x,ct_y} = [];
                Parameters{ct_x,ct_y}      = [];
            end
            
        end
    end
    
    % *********************************************************************
    % Save scaling factors
    % *********************************************************************
    file_save = [SAT_SST_func_IO('data',P),'/Valiation_20_year_statistics_',P.date,'.mat'];
    save(file_save,'R','E','-v7.3');

    file_save = [SAT_SST_func_IO('data',P),'/Valiation_20_year_data_',P.date,'.mat'];
    save(file_save,'Validation_data','Parameters','-v7.3');
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Out_training,Out_validation] = prepare_grid_data(Out,ct_x,ct_y,P)

    l = Out.lon > (ct_x*10-10)  & Out.lon <= (ct_x*10) & ...
        Out.lat > (ct_y*10-100) & Out.lat <= (ct_y*10-90);

    if nnz(l) > 0
        disp([ct_x ct_y])
        Out_grid = pre_process_data(Out,l,P);
        
        % split training and validation data ------------------------------
        if ~isempty(Out_grid)
            Out_training{1}                         = Out_grid;
            Out_training{1}.AT(1:240)               = [];
            Out_training{1}.SST(1:240)              = [];
            Out_training{1}.AT_itp(1:240)           = [];
            Out_training{1}.SST_itp(1:240)          = [];
            Out_training{1}.l_use_loss(1:240)       = [];
            Out_validation{1}                       = Out_grid;
            Out_validation{1}.AT(241:end)           = [];
            Out_validation{1}.SST(241:end)          = [];
            Out_validation{1}.AT_itp(241:end)       = [];
            Out_validation{1}.SST_itp(241:end)      = [];
            Out_validation{1}.l_use_loss(241:end)   = [];
            
            Out_training{2}                         = Out_grid;
            Out_training{2}.AT(end-239:end)         = [];
            Out_training{2}.SST(end-239:end)        = [];
            Out_training{2}.AT_itp(end-239:end)     = [];
            Out_training{2}.SST_itp(end-239:end)    = [];
            Out_training{2}.l_use_loss(end-239:end) = [];
            Out_validation{2}                       = Out_grid;
            Out_validation{2}.AT(1:end-240)         = [];
            Out_validation{2}.SST(1:end-240)        = [];
            Out_validation{2}.AT_itp(1:end-240)     = [];
            Out_validation{2}.SST_itp(1:end-240)    = [];
            Out_validation{2}.l_use_loss(1:end-240) = [];
        else
            Out_training   = [];
            Out_validation = [];
        end
    else
        Out_training   = [];
        Out_validation = [];
    end
end

function Out_grid = pre_process_data(Out,l,P)

    Out_grid.N   = nnz(l);

    Out_grid.lf  = nanmean(Out.land_fraction(l));

    AT  = CDC_pnt2glbmean(Out.raw_SAT_anm(l,:,:),Out.lon(l),Out.lat(l));
    SST = CDC_pnt2glbmean(Out.raw_SST_anm(l,:,:),Out.lon(l),Out.lat(l));

    yrs = P.yr_learn - 1849;
    AT  = AT(:,yrs);  AT  = AT(:);
    SST = SST(:,yrs); SST = SST(:);

    if nnz(~isnan(AT + SST)) > 240

        % Subset years where both AT and SST has values ---------------
        temp = reshape(AT+SST,12,numel(P.yr_learn));
        l_yr = any(~isnan(temp),1);
        l_yr(find(l_yr,1):find(l_yr,1,'last')) = 1;
        yrs  = P.yr_learn;
        Out_grid.yrs  = yrs(l_yr);

        % Linear interpolate if missing values ------------------------
        l_yr = repmat(l_yr,12,1);
        Out_grid.AT   = AT(l_yr);
        Out_grid.SST  = SST(l_yr);

        Out_grid.AT_itp  = CDC_interp1(Out_grid.AT,[],1);
        Out_grid.SST_itp = CDC_interp1(Out_grid.SST,isnan(Out_grid.AT),1);

        Out_grid.l_use_loss  = ~isnan(Out_grid.AT) & ~isnan(Out_grid.SST);
    else
        Out_grid = [];
    end
end

function Out = read_data(P)

    % *********************************************************************
    % Load data
    % *********************************************************************
    file_load = [SAT_SST_func_IO('data',P),'/Raw_SAT_SST_',P.date,'.mat'];
    load(file_load,'raw_SST_anm','lon','lat','land_fraction');

    file_load = [SAT_SST_func_IO('data',P),'/Raw_SAT_anm_',P.date,'.mat'];
    load(file_load,'raw_SAT_anm')

    % *********************************************************************
    % Subset long stations
    % *********************************************************************
    % GLB_mod_subset_long_stations;

    raw_SST_anm     = raw_SST_anm - nanmean(raw_SST_anm(:,:,[1982:2014]-1849),3);
    raw_SAT_anm     = raw_SAT_anm - nanmean(raw_SAT_anm(:,:,[1982:2014]-1849),3);

    Out.raw_SST_anm = raw_SST_anm;
    Out.raw_SAT_anm = raw_SAT_anm;
    Out.lon         = lon;
    Out.lat         = lat;
    Out.land_fraction = land_fraction;
    
end
